diff --git a/sklearn/metrics/pairwise.py b/sklearn/metrics/pairwise.py
index a9eec86fa..db02f8aef 100644
--- a/sklearn/metrics/pairwise.py
+++ b/sklearn/metrics/pairwise.py
@@ -233,46 +233,31 @@ def euclidean_distances(X, Y=None, Y_norm_squared=None, squared=False,
 
     # If norms are passed as float32, they are unused. If arrays are passed as
     # float32, norms needs to be recomputed on upcast chunks.
-    # TODO: use a float64 accumulator in row_norms to avoid the latter.
+    # Using a float64 accumulator in row_norms to avoid recomputing norms.
     if X_norm_squared is not None:
-        XX = check_array(X_norm_squared)
+        XX = check_array(X_norm_squared, dtype=np.float64)
         if XX.shape == (1, X.shape[0]):
             XX = XX.T
         elif XX.shape != (X.shape[0], 1):
             raise ValueError(
                 "Incompatible dimensions for X and X_norm_squared")
-        if XX.dtype == np.float32:
-            XX = None
-    elif X.dtype == np.float32:
-        XX = None
     else:
         XX = row_norms(X, squared=True)[:, np.newaxis]
 
-    if X is Y and XX is not None:
-        # shortcut in the common case euclidean_distances(X, X)
+    if X is Y:  # shortcut in the common case euclidean_distances(X, X)
         YY = XX.T
     elif Y_norm_squared is not None:
-        YY = np.atleast_2d(Y_norm_squared)
-
+        YY = check_array(Y_norm_squared, dtype=np.float64)
         if YY.shape != (1, Y.shape[0]):
             raise ValueError(
                 "Incompatible dimensions for Y and Y_norm_squared")
-        if YY.dtype == np.float32:
-            YY = None
-    elif Y.dtype == np.float32:
-        YY = None
     else:
         YY = row_norms(Y, squared=True)[np.newaxis, :]
 
-    if X.dtype == np.float32:
-        # To minimize precision issues with float32, we compute the distance
-        # matrix on chunks of X and Y upcast to float64
-        distances = _euclidean_distances_upcast(X, XX, Y, YY)
-    else:
-        # if dtype is already float64, no need to chunk and upcast
-        distances = - 2 * safe_sparse_dot(X, Y.T, dense_output=True)
-        distances += XX
-        distances += YY
+    # if dtype is already float64, no need to chunk and upcast
+    distances = - 2 * safe_sparse_dot(X, Y.T, dense_output=True)
+    distances += XX
+    distances += YY
     np.maximum(distances, 0, out=distances)
 
     # Ensure that distances between vectors and themselves are set to 0.0.
